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The growing interest in carbon-based spintronics has stimulated a number of recent theoretical 
studies on the RKKY interaction in graphene, based on which the energetically favourable alignment 
between magnetic moments embedded in this material can be calculated. The general consensus is 
that the strength of the RKKY interaction in undoped graphene decays as 1/D 3 or faster, where 
D is the separation between magnetic moments. Such an unusually fast decay for a 2-dimensional 
system suggests that the RKKY interaction may be too short ranged to be experimentally observed 
in graphene. Here we show in a mathematically transparent form that a far more long ranged 
interaction arises when the magnetic moments are taken out of their equilibrium positions and 
set in motion. We not only show that this dynamic version of the RKKY interaction in graphene 
decays far more slowly but also propose how it can be observed with currently available experimental 
methods. 

PACS numbers: 



I. INTRODUCTION 

The atomically thin sheet of carbon known as graphene 
has been attracting the interest of the wider scientific 
community due to its enormous potential for applications 
in fields as diverse as photonics, sensor technology and 
spintronics, to name but a few-ir— . Spintronics is a par- 
ticularly promising field for graphene application due to 
the weak spin-orbit and hyperfine interactions, which in 
other materials act as significant sources of spin relax- 
ation and decoherence^— 

One recurrent topic in the field of spintronics is the 
mechanism by which localized magnetic moments em- 
bedded into nanoscale systems are able to interact with 
each other even when far apart. An indirect exchange in- 
teraction mediated by the conduction electrons of a host 
medium manifests itself as an energy difference between 
different alignments of the localized moments, leading 
to energetically favourable configurations. Such an in- 
teraction is usually calculated within the Ruderman- 
Kittel-Kasuya-Yosida (RKKY) approximation^— and 
indeed the interaction itself often takes this moniker. 
The RKKY interaction in graphene has been intensively 
studied^— and there is a general consensus that its 
strength decays faster than in other 2-dimensional ma- 
terials. In fact, the RKKY interaction between magnetic 
moments a distance D apart decays as 1/D 3 or faster 
in undoped graphene. Such a fast decay rate indicates 
that the interaction is rather short ranged and possibly 
not as easily observable as, for instance, in the case of 
carbon nanotubes^, for which the RKKY interaction is 
predicted^— to decay as 1/D. In the case of doped 
graphene^ the RKKY interaction decays as 1/D 2 , as 
commonly expected for 2-dimensional materials, but it 
still remains to be experimentally seen. 

Rather than ruling out the RKKY in graphene as too 
short ranged an interaction to be seen in practice, here 



we show that its range is considerably augmented if the 
magnetic moments are set in motion. The notion of 
a dynamic RKKY has been proposed by Simanek and 
Hcinrich^ who generalized the concept of the electron- 
mediated interaction between magnetic moments to the 
case of when these moments are not in static equilibrium. 
Moreover, an expanded range of the magnetic interac- 
tion has been reported in multilayered materials when 
their magnetization varies with time^i, although no di- 
rect link with the RKKY interaction was established. In 
this manuscript we demonstrate that the magnetic inter- 
action range in graphene can be expanded in a similar 
fashion. Mathematically transparent expressions for the 
dynamic spin susceptibility of this material are derived, 
based on which we identify a striking modification in the 
decay rate of the RKKY interaction when the magnetic 
moments are allowed to precess. 



II. MODEL 

Let us start by recalling that the conventional RKKY 
interaction between two magnetic objects separated by 
a distance D results from the energy involved in placing 
one of the objects in the presence of the electronic spin- 
polarization induced by the other. When minimizing this 
energy, the most favourable alignment of the magnetiza- 
tions is found. It is thus not surprising that the strength 
of the RKKY decays with D in exactly the same way 
as the induced spin polarization due to a single object. 
In other words, when determining the decay rate of the 
RKKY interaction, it is sufficient to consider how the in- 
duced spin polarization caused by a single magnetic ob- 
ject is spread in space. With that in mind, we shall first 
consider a single magnetic object on a graphene sheet. 
For the sake of simplicity, let us assume that the object 
is represented by a single substitutional magnetic atom. 
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This assumption can be relaxed to consider more gen- 
eral magnetic objects without significantly affecting the 
conclusions obtained. 

Instead of studying the effect of this magnetic atom 
in a static situation, we impose a perturbing time- 
dependent transverse magnetic field of the form 



ho [cos(wt) x — sin(u;i) y] 



(1) 



that sets the moment in precession with a frequency u>. 
Experimentally this precession can be achieved in a num- 
ber of different ways^~— . The spin polarization induced 
by the precessing moment becomes time dependent so 
that a second magnetic object would couple directly with 
the oscillatory polarization and precess in exactly the 
same way. This polarization is expressed by the trans- 
verse spin susceptibility x, which reflects how the spin 
degrees of freedom of a system respond to a magnetic 
excitation. 

To calculate the spin susceptibility one needs the 
Hamiltonian describing the electronic structure of the un- 
perturbed system, which we assume is given by 

H = ^2 7i? 4a +y^( £ 0"-0,g + ~ flQ.a WQ,ff) + H Z ■ 
<i,j>,<J a 

Here, ^ represents the electron hopping between near- 
est neighbor sites i and j, c| ff (cj CT ) creates (annihilates) 
an electron with spin a in site i, eo is the atomic en- 
ergy level of the magnetic site chosen to be located at 
the origin (site 0), fiQ a = Cq a co, a is the corresponding 
electronic occupation number operator, and U represents 
an effective on-site interaction between electrons on the 
magnetic site, which is neglected elsewhere. Finally, Hz 



plays the role of a local Zeeman interaction that defines 
the z-axis as the equilibrium direction of the magneti- 
zation. This unperturbed Hamiltonian is equivalent to 
that of the Anderson model used to describe localised 
magnetic impurity states in metals^. The interaction 
between the oscillatory magnetic field in Eq. ([lj and the 
magnetization at site 0, So, is accounted for by an inter- 
action term Hi nt = g/ishj^ ■ So- Our results are not crit- 
ically dependent on the electronic structure parameters 
of the Hamiltonian, indicating that our findings are valid 
for a variety of magnetic objects. In fact, as far as the 
decay rate is concerned, it is fully determined by the host 
material, in this case the graphene sheet. Finally, spin- 
orbit coupling is neglected due to the long spin-diffusion 
length in graphene. 

The time-dependent transverse spin susceptibility is 
defined as xy(t) = -|6(t)([-S+(t), ^"(0)]), where e(x) 

is the Heaviside step function, and 5+ and S~ are the 
spin raising and lowering operators at site m, respec- 
tively. The indices j and m refer to the locations where 
the field is applied and where the response is measured, 
respectively. In our case, a precession of the magnetic 
moment is induced at the magnetic site 0, and we wish 
to observe the spin disturbance a distance D apart at 
an arbitrary site m. This response is fully described by 
Xm,o{t)- Within the random phase approximation, this 
susceptibility may be calculated in the frequency domain 
as 



(3) 



where \° is the Hartree-Fock susceptibility, whose matrix 
elements are given by 



J 



■ t /*+oo 

xlM = ^ / dE'f{E>) { [Gl,{E') GjJ(E')] G\.{E' + M + [cf^E') - G^{E')\ Gj}(E> m} ■ (4) 
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Here, GJj(E') and GjJ(E') represent the time Fourier 
transforms of the retarded and advanced single-electron 
Green Functions (GF), respectively, for an electron of en- 
ergy E' with spin a between sites I and j, and f(E') is 
the Fermi function. Concerning the physical meaning of 
the quantity Xm,o(w), its absolute value is proportional 
to the amplitude of the spin precession at site m as a 
response to a precession induced at site 0. Therefore, to 
assess how a second magnetic object inserted at site m 
will couple with that located at the origin, the relevant 
matrix element of the susceptibility is Xm,o- To simplify 
the analytic calculation of the interaction behaviour we 
consider the case of a single impurity. However, the ap- 
proach can be generalised to the cases of magnetic clus- 
ters or multiple impurities by introducing a sum over im- 



purity sites in the second term of the Hamiltonian given 
by Eq. ©. Such an approach is used in numerical cal- 
culations involving more than one impurity later in this 
work. 



From Eqs. ([3]) and (@| it is clear that the position de- 
pendence of the spin susceptibility is entirely contained 
within the corresponding Hartree-Fock term \m o- Fur- 
thermore, the dependence of this term on the separation 
D between the sites and m can be calculated more 
easily by replacing the matrix elements of G a (E') that 
appear in Eq.Q with the bulk GF of pristine graphene, 
which we calculate within the tight-binding formalism. 
In this case Xm.o( w ) i s written as \m o( w ) = I\ + h + h, 
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when 
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III. INTERACTION RANGE 



h = 7T / dE'f(E')g 0}m (E')g mt0 (E' + huj) 
h = -5- / dE'f(E')g- fi (E')g^ m (E' - M 

^ J-oo 

i f +co 

h = dE ' if( E + M - f( E ')] 5o,m( E 'Wm,o( E ' 

(5) 

Here g (g~) is the retarded (advanced) GF of pristine 
graphene. Note that Ii (J 2 ) involves the convolution of 
two retarded (advanced) GF whereas I3 involves one of 
each type. 

We can now make use of the simplicity of the electronic 
structure of graphene and write the off-diagonal matrix 
elements of the pristine retarded GF of graphene as 



g , m (E) = g m , ( E ) = 



A(E) e* 



Q(E)D 



D 



(0) 



This expression for the off-diagonal GF is valid for mod- 
erately large values of D.— Furthermore, it is valid for 
energies not only around the Fermi level where the dis- 
persion relation is linear but also across the entire energy 
band. The precise forms of the functions A(E) and Q(E) 
depend on the direction of the vector joining the sites 
and m on the graphene sheet but what is obvious from 
the expression is the simple functional form of 5o,m as 
a function of the separation D. For the case where the 
vector joining the sites and m is along the armchair di- 
rection and both sites are located on the same sublattice, 
we can write 



Q(z) = iskr 1 



and 




(z 2 + 3 7 2 ) V7 2 - z 2 ' 



(7) 



(8) 



where 7 is the nearest neighbour hopping of graphene. 
In the next section we investigate how the amplitude of 
the dynamic susceptibility varies as the separation is in- 
creased along this direction. Changing the direction in- 
vestigated may introduce additional oscillatory features, 
similar to those previously noted in the static case, but 
will not affect the decay rate of the interaction, which 
is our primary concern in this work. Allowing impuri- 
ties to locate on different sublattices leads to qualitative 
changes in the static interaction where the sign of the 
coupling indicates preferred magnetic alignments. How- 
ever, once again, the decay rate is not altered. This will 
be seen in Section [IV! where multiple impurities located 
on both sublattices and separated in random directions 
are considered numerically to reproduce a more likely ex- 
perimental configuration. 



We can now investigate the decay rate of the RKKY 
interaction by studying the three integrals in Eq. ([5]) sep- 
arately. By symmetry, I\ and I2 have similar functional 
forms and therefore identical contributions. For the sake 
of conciseness, let us focus on the integral I\ . Combining 
fj&fs.- © an d @ we may write 



+ °° B(E' u})e l(Q(E ^ +Q( - E ' +hu ^ D 

dE' [ ' J „,„, - „ , (9) 



D(l 



where E F is the Fermi energy, B(E' , uj) = A(E')xA(E' + 
hio), (3 = jr-jt, T being the temperature and ks the 
Boltzmann constant. The integral in Eq © can be solved 
by employing a semi-circular contour in the upper-half 
of the complex energy plane to reduce the integral to a 
sum over Matsubara frequencies. Expanding the func- 
tions B(E',ui) and Q(E) around Ep and taking the low 
temperature limit, T — > 0, we find that 1\ is given by22 



t = 1 i\2Q(E F ) + Q'(E F )hu]D \p (~iy +1 B e (E F ,0) 1 

1 2tt 4^ [2iQ'{Ef)Y +1 D l + 2 



(10) 



where I is an integer and B l {Ep,Q) is the corresponding 
l th -order derivative of the function B(E' , 0) evaluated at 
Ef , which results from its Taylor expansion. 

The Z?-dependence of I\ is now evident. In the asymp- 
totic limi t 22 ^ 31 of large separations it is determined by the 
leading term in Eq. (fTT)|) . namely t = 0, suggesting that, 
in general, I\ ~ D~ 2 . However, at Ep = 0, the coefficient 
£?(°) vanishes and the decay rate is in fact determined by 
the first surviving term, I = 1, resulting in l\ ~ D~ 3 for 
undoped graphene. When Ep ^ 0, B^ does not vanish, 
and in this case I\ does decay as D~ 2 . Notice that the 
frequency-dependence is contained only in the argument 
of the exponential in Eq. ([I0[) . meaning that even when 
cj = the decay rate of I\ remains unaltered. This is 
precisely what is found for the static RKKY and agrees 
with the aforementioned feature that this interaction de- 
cays as D~ 3 (D~ 2 ) for undoped (doped) graphene. 

Regarding the -Z3 integral, for ksT — > 0, the step-like 
Fermi functions can be used to rewrite this term as 



h = 



dE'gi m (E')g , m (E' + hj) . (ii) 



The pristine Green function <7o,m given by Eq.© can 
be more concisely expressed for the small integration 



range around Ep if one rewrites Q(z) 



and A(z) 



respectively. In this case, the distance de- 



pendence can be taken outside the integral and we find 

h 



3vrL>|7| 



dE' yf-E'{E' + hu) 



EF — fibJ 



(12) 
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the reported RKKY interaction that goes as D~ 3 (D~ 2 ) 
when Ep = (Ep ^ 0). Attention is drawn to the 1/D 
asymptotic decay that appears in all dynamic cases and 
which indicates that the RKKY should indeed be more 
long ranged when the magnetic moments are in motion. 
This behavior is further confirmed by the respective in- 
sets which correspond to calculations of the spin suscep- 
tibility for the case of two magnetic impurities a distance 
D apart. For these calculations the susceptibility in Eq. 
(U) uses numerical Green functions whose corresponding 
Hamiltonian includes the additonal terms necessary to 
describe a second magnetic impurity. In this case the 
spin disturbance at the second moment induced by the 
precession of the first is calculated. Note that the respec- 
tive decay rates in this case are in excellent agreement 
with the predictions made for a single impurity, confirm- 
ing the analytic treatment of the decay rate behaviour 
outlined above. 



IV. EXPERIMENTAL SIGNATURES 



FIG. 1: Log-log plots of the numerically calculated spin sus- 
ceptibility (in arbitrary units) probed at a distance D from 
a single magnetic impurity. Top panels correspond to the 
case of static susceptibility whereas bottom panels are for fi- 
nite frequencies. Left (Right) panels are for undoped (doped) 
cases. (Red) Dots are the numerically calculated values and 
(black) lines display the predicted power laws also highlighted 
on each of the main panels. Insets display the corresponding 
calculations for two impurities where exactly the same power 
law behavior is found. 



which in the undoped (Ep = 0) case reduces to 



h = - 



2AttD\j[ 



(13) 



The separation dependence coming from /3 is clearly 
more long ranged than that from I\ and I2 since Eq. (fT2)) 
decays as D . Furthermore, unlike the expression for l\ 
in Eq. (|10j) which remains finite even in the static limit, 
the contribution from I3 vanishes when u = 0. There- 
fore, the RKKY interaction between localized magnetic 
moments in graphene becomes more long ranged once 
the moments are set in motion with a finite excitation 
frequency. 

Figure 1 confirms the analytically predicted decay 
rates with a fully numerical evaluation of the spin sus- 
ceptibility using Eqs. <j3j> and (j4}. The top panels depict 
IXm,o| for a single magnetic impurity located at site in a 
purely static configuration, i.e. lu = 0, whereas the bot- 
tom panels show the same result for the case in which 
the moment is precessing at the resonance frequency. 
Left (right) panels correspond to undoped (doped) cases. 
(Red) Dots are the numerically evaluated results and the 
(black) lines follow the power laws shown in the Figure. 
In the static case the spin polarization induced by the 
magnetic object decays as expected, at the same rate as 



To date it has been very difficult to probe the RKKY 
interaction in graphene experimentally. It is understand- 
able that with a decay rate as fast as D~ 3 it is difficult 
to probe the interaction for any reasonable separation. 
The presence of magnetism in disordered graphene sys- 
tems may indicate the presence of an exchange coupling 
between magnetic moments formed around defects. Nu- 
clear magnetic resonance experiments reveal that these 
defects have indeed magnetic moments, since they couple 
to implanted Fe atoms^. However, whether or not these 
moments couple with each other, or with the graphene 
lattice, to form a ferromagnetic state is a controversial 
subject and many of the results in this area have proved 
difficult to reproduce 3 . However, our results suggest 
that magnetic moments should be able to feel their mu- 
tual presence at greater separations once they are set 
in motion. To address the issue of how this dynamic 
RKKY interaction can be probed, we turn our attention 
to a slightly different feature that can also be extracted 
from the spin susceptibility, namely the lifetime of mag- 
netic excitations, something that can be probed using 
the method of inelastic scanning tunneling spectroscopy 
(ISTS)4irii Rather than looking at the off-diagonal el- 
ement of the susceptibility for a given frequency, the di- 
agonal matrix element Xm.m plotted as a function of the 
frequency uj tells us how strongly the system responds 
to a time-dependent magnetic excitation. In particu- 
lar, peaks in the w-dependent susceptibility reflect the 
existence of resonance frequencies whereas their inverse 
widths characterize the respective lifetime of the spin ex- 
citations. Figure 2A shows xo.o a $ a function of the exci- 
tation frequency u for a single magnetic impurity at the 
origin (site 0), where a single peak is clearly identifiable 
with a linewidth W\ . 

To examine how the dynamic RKKY interaction be- 
tween moments affects the lifetime of the spin excita- 
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FIG. 2: (A) Diagonal element of the spin susceptibility xo,o 
(in arbitrary units) obtained for a single magnetic impurity 
located at site plotted as a function of the excitation fre- 
quency ui. A single peak of width W\ is highlighted. (B) 
Fractional deviations of the linewidth ((W — Wi)/Wi) plot- 
ted as a function of the relevant separation (in units of the 
lattice parameter of graphene). The solid (blue) line corre- 
sponds to the case of two impurities a distance D apart. The 
gray area delimits the region spanned by the two- impurity re- 
sult if all possible directions are included. (C) The scattered 
(red) dots correspond to the linewidths W of the susceptibil- 
ity Xm.m for different sites m in the case of a disordered array 
of magnetic impurities randomly distributed on graphene, as 
shown schematically in the inset, plotted as a function of the 
corresponding nearest neighbor separation, Dmn, i.e. the dis- 
tance to the impurity nearest m in the configuration. 



tions, a second moment is introduced into the system as 
before and its distance, D, from the first moment varied 
along the armchair direction. We have seen clearly in the 
previous section how the interaction between precessing 
moments manifests itself in the off-diagonal term of the 
spin susceptibility. It is now worth investigating how the 
diagonal term is affected since the resultant variations 
in the spin excitation spectrum may provide experimen- 
tally detectable signatures of a dynamic RKKY interac- 
tion. We expect that when the moments are moved far 
enough apart they become essentially independent and 
the excitation linewidth should approach that of the sin- 
gle impurity case. In fact, we can use Eq. ([3]) to expand 
the diagonal term of the susceptibility for two magnetic 



(2) 

impurities at sites and to, x q, as 



(2) 



Y (l) 

AO.O 



{i + x%u) (1 



Am,m 



(14) 



where x'oq is the single-impurity susceptibility. The first- 
order correction to the diagonal matrix element of the 
single-impurity susceptibility when a second moment is 
introduced at site m contains the product Xo m Xm o an< ^ 
should therefore decay as D~ 2 , with higher order terms 
in the expansion decaying even faster. This analysis sug- 
gests that fluctuations in the magnetic excitation spec- 
trum induced by the introduction of an additional mo- 
ment should decay as D~ 2 towards the single-impurity 
result. 

Our numerical simulations find that adding a second 
impurity to the system maintains the peaked features 
in the spin susceptibility found for the single impurity 
case in Figure 2A, but that now the linewidths become 
Independent. This can be viewed in Figure 2B where 
AW/Wi represents the relative deviation of the linewidth 
W for the two impurity case from that of the single- 
impurity linewidth, W\ . The solid (blue) line shows that 
the linewidth oscillates around the value for the one- 
impurity case, with an amplitude that decays asymp- 
totically towards it as the separation between the two 
moments is increased. This decay rate is found to be 
D~ 2 , as predicted above and indicated by the shaded 
area enveloping the two-impurity curve. It is worth not- 
ing that a mere phase shift in the solid line occurs if one 
changes the direction along which the two impurities are 
separated. If all possible directions were considered, the 
resulting data points would generate the shaded area in 
Figures 2B and 2C. 

Rather than having one or two isolated impurities, a 
more realistic scenario is the case depicted in the inset 
of Figure 2C which shows several magnetic atoms ran- 
domly located across the graphene sheet. Many disor- 
dered configurations of impurities are generated where 
the impurities are allowed to locate in arbitrary direc- 
tions from each other and occupy either sublattice. For 
each magnetic impurity we calculate the linewidth of the 
diagonal spin susceptibility and also record the separa- 
tion from the nearest other impurity in that configura- 
tion. These linewidths are plotted as relative deviations 
as before, this time as a function of nearest-neighbour dis- 
tance, Datat, and appear in the main panel of Figure 2C 
as the solid (red) scattered dots. Notice that the spread 
of calculated linewidths for the disordered configurations 
agrees well with the shaded area that enveloped the 
two-impurity linewidth oscillation. Therefore, the dy- 
namic interaction between the nearest-neighbour impuri- 
ties tends to dominate and determine the overall lifetime 
of the spin excitation in these disordered configurations. 
Consequently, the fluctuations in the linewidth measure- 
ments reflect the variations of the dynamic RKKY in- 
teraction and the spread of measurements for a given 
distance, i.e. the shaded area, decays as D~ 2 . Bearing 
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in mind that the nearest-neighbour separation is directly 
related to the concentration of impurities p ~ Djj 2 N , we 
predict that the standard deviation (aw) of the mea- 
surements of the spin-excitation lifetimes will scale with 
the concentration of magnetic impurities as a\y ~ p. 
This is an unmistakable signature of the RKKY inter- 
action in its dynamic form. Note that the inset of Figure 
2C matches the experimental setup of Reference 4 ^ which 
measured the spin excitation lifetimes of randomly dis- 
persed Fe atoms on Cu surfaces with ISTS. It should be 
equally possible to probe the dynamic RKKY interac- 
tion in graphene using similar techniques since C atoms 
have substantially smaller spin-orbit coupling than Cu. 
We note that our studies have assumed low concentra- 
tions of magnetic impurities and also atomically precise 
graphene samples, which may be difficult to achieve ex- 
perimentally. A recent study suggests that, for the static 
RKKY, strong on-site disorder can induce an exponen- 
tial suppression of the coupling at large distances^. We 
expect that similar features may occur in the dynamic 
RKKY under the same kind of disorder regimes. 

V. CONCLUSION 

In summary, we have shown that the RKKY inter- 
action between localized magnetic moments embedded 
in a non-magnetic material, which is predicted to de- 
cay rather fast in graphene, may become substantially 
more long ranged once the magnetic moments are taken 



out of equilibrium and set to precess. This has been il- 
lustrated using analytic arguments and fully numerical 
calculations. We argue that the difficulty of experimen- 
tally probing the RKKY interaction in graphene-based 
materials may be overcome by simply exciting the mag- 
netization of magnetic objects in contact to graphene, 
something that is currently achievable with ISTS mea- 
surements. Finally, because our conclusions on the in- 
creased range of the RKKY interaction are equally valid 
for a variety of other magnetic objects, such as adatoms, 
nanoparticles and vacancy-induced magnetic moments, 
to name but a few, this may open the road to new ways 
of transporting magnetic information in a material that 
carries very weak spin-orbit coupling and which thereby 
dissipates very little spin current. 
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